function y_sub = expsub(y_data, exp_order);

% EXPSUB - subtract an exponential fit of the data
%
%    YS = EXPSUB(Y, N) subtracts an N-exponential fit of the
%    rows in Y from Y and stores the result in YS.

% By:   S.C. Molitor (smolitor@med.unc.edu)
% Date: September 29, 2000

% initialize output

y_sub = [];

% validate input arguments

if (nargin ~= 2)
   msgbox('Invalid number of arguments', 'EXPSUB Error', 'warn');
   return
elseif (~isnumeric(y_data) | isempty(y_data))
   msgbox('Y must be a numeric array', 'EXPSUB Error', 'warn');
   return
elseif (~isscalar(exp_order))
   msgbox('N must be a scalar', 'EXPSUB Error', 'warn');
   return
elseif ((exp_order < 1))
   y_sub = y_data;
   return
end

% initialize exponential fit parameters

x_data = [0 : size(y_data, 2) - 1];
y_steady = mean(y_data(:, end - floor(size(y_data, 2)/10) : end), 2);
y_amp = y_data(:, 1) - y_steady;
for i = 1 : size(y_data, 1)
   y_tau(i) = sum(y_data(i, :) - y_steady(i)) / y_amp(i);
end

% account for multiple exponentials

y_par = [y_steady y_amp abs(y_tau)'];
for i = 2 : exp_order
   y_amp = y_amp/5;
   y_tau = 5*y_tau;
   y_par(end - 1, :) = y_par(end - 1, :) - y_amp;
   y_par = [y_par y_amp y_tau];
end

% loop to fit exponentials
% setup & update waitbar

y_sub = y_data;
%h_wait = waitbar(0, '0% completed...');
%h_axes = findobj(h_wait, 'Type', 'axes');
%set(h_wait, 'Name', 'Exponential Subtraction Progress');
for i = 1 : size(y_data, 1)
   [p, chisq, niter, y_fit] = mrqfit_scm('exponential', y_par(i, :), x_data, y_sub(i, :));
   y_sub(i, :) = y_sub(i, :) - y_fit;
   %waitbar(i/size(y_data, 1));
   %set(get(h_axes(1), 'Title'), 'String', sprintf('%d%% completed...', round(100*i/size(y_data, 1))));
   %drawnow;
end
%delete(h_wait);
return
